The Thermodynamic Stability of Two Dimensional Crystals with an Extended 

Coupling Scheme 
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We calculate mean square deviations for crystals in one and two dimensions. For the two di- 
mensional lattices, we consider several distinct geometries (i.e. square, triangular, and honeycomb), 
and we find the same essential phenomena for each lattice structure. We investigate the stability of 
long-range crystalline order for a variety of coupling schemes, including short-range exponentially 
decaying inter-atomic potentials and long-range interactions with a power law dependence r~ a . For 
the latter in the ID case, we find a critical value a™ = 1.615 ± 0.005 for the power law decay 
exponent below which crystalline order is intact, and above which thermal fluctuations destroy 
long-range order when T > 0. The corresponding critical value for two dimensional lattices with 



displacements confined to the plane is cq. 



3.15 ± 0.025. If motion perpendicular to the crystal 



plane is permitted, thermally induced distortions diverge rapidly (i.e. linearly) in dual layer systems 
with local stiffness provided by an extended coupling scheme, even if the interaction is long ranged, 
decaying as a power law in the separation between lattice sites. 



PACS numbers: 62.25.Jk, 62.23.Kn, 63.22.Np 



I. INTRODUCTION AND THEORETICAL 
TECHNIQUES 

Crystals are regular periodic physical systems with the 
atomic constituents organized in periodic arrays, where 
the periodicity is a characteristic of the crystal in its equi- 
librium configuration and a manifestation of long-range 
positional order. However, the effect of thermal fluctu- 
ations must be taken into consideration at finite tem- 
peratures where thermally excited lattice vibrations may 
degrade or destroy crystalline order. An early theoreti- 
cal treatment developed by Lindemann 1 examined the ef- 
fect of lattice vibrations in a framework which neglected 
the correlations of atomic motions, but which nonetheless 
provides a reasonable description (on an order of magni- 
tude basis) for the melting of three dimensional solids. 

A salient component of the Lindemann analysis is the 
Lindemann criterion where the melting of a solid is con- 
sidered to have taken place when mean square deviations 
from equilibrium exceed a tenth of a lattice constant. X- 
ray diffraction data, which may provide a direct measure 
of long-range order in a crystal lattice (and hence a means 
to determine temperatures where crystalline order is lost) 
finds reasonable agreement^ with the Lindemann crite- 
rion. The accord is manifest in the finding that Bragg 
peaks corresponding to broken translational symmetry 
vanish when mean square fluctuations from equilibrium 
(also determined from an analysis of X-ray diffraction 
data) are in the vicinity of a tenth of a lattice constant, 
as specified in the Lindemann result. 

A factor of significance for the effectiveness of the Lin- 
demann criterion is the tendency for atoms in three di- 
mensional crystals to have a large number of neighbors 
[e.g. a dozen nearest neighbors in the face of face centered 
cubic (fee) lattices] . Hence, mean field treatments in the 
spirit of Weiss molecular mean field theory are more likely 
to provide a reasonable theoretical description since sta- 



tistical fluctuations tend to suppressed somewhat by av- 
eraging when a relatively large number of neighbors are 
present. 

On the other hand, it should be understood that apart 
from the number of nearest neighbors, dimensionality is 
a very important parameter which may affect the ther- 
modynamic characteristics and integrity of a crystal lat- 
tice to a large degree. Nano-scale engineering often takes 
place in systems of low dimensionality such as carbon- 
nanotubes where the length may exceed the width by 
several orders of magnitude; nanotubes tend to be re- 
garded as one dimensional systems. Graphene sheets, 
covalently bonded single layer honeycomb lattices of car- 
bon atoms, may be considered genuine monolayers, and 
hence possess strongly two dimensional character. The 
thermodynamic stability of a system is strongly depen- 
dent on its dimensionality with statistical fluctuations 
becoming more important for two dimensional systems, 
and very important for essentially ID structures such as 
nanotubes. 

An important theoretical result known as the Mermin- 
Wagner theorem^ predicts that as the bulk limit is ap- 
proached, thermal fluctuations destroy long-range crys- 
talline order in the context of ID lattices. However, this 
result does not preclude the stabilization of positional 
order for T > if the interaction between atomic mem- 
bers is long-ranged (e.g. decaying as a power law in the 
separation between positions in the crystal lattice). In 
this work, it is our program to examine conditions which 
preserve long range order in low dimensional systems at 
finite temperatures. 

In one dimension, the deleterious effect of thermal fluc- 
tuations is felt most severely, and ultimately only short- 
range order exists if the interatomic interaction is finite 
in range. In three dimensions, long-range positional or- 
der is intact for finite temperatures below the melting 
temperature T m . Two dimensional solids are often re- 
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garded as an intermediate case where thermal fluctua- 
tions are strong enough to destroy long-range order as 
the size of the system is increased, but only in a very 
gradual manner. Although crystalline order in 2D sys- 
tems does not survive in the thermodynamic limit if the 
interaction is confined to nearest neighbors or is other- 
wise finite in range, a long-range coupling with power law 
decay may stabilize long-range positional order. In fact, 
even for one dimensional solids, we find a critical decay 
exponent a]P — 1.615 ± 0.005 below which crystalline 
order remains stable for T > 0, whereas long-range order 
is only gradually lost in the bulk for power law decays 
where a > a c . Similarly, the corresponding exponent in 
2D is a^P ~ 3.15 ± 0.025, where long range crystalline 
order is preserved for a < a^ D , whereas thermal fluctu- 
ations destroy positional order if a > a^°. Within the 
bounds of numerical error, we obtain the same value for 
the threshold exponent a): D for distinct lattice geometries 
including square lattices, triangular lattices, and honey- 
comb lattices. 

We examine various types of coupling schemes, includ- 
ing very short-ranged interactions where atoms interact 
with only a few nearest neighbors and perhaps also next- 
nearest neighbors. We also consider extended schemes 
where there is a finite coupling to all neighbors, but where 
the interaction is still short-ranged, with a rapid decay 
profile, such as that of an inter-atomic potential with an 
exponential dependence V cx e _7r where 7 -1 is the fi- 
nite length scale corresponding to the coupling scheme. 
Finally, we also consider a long-ranged algebraically de- 
caying coupling of the form V(r) = r~ a where a may 
assume different values, though for the energy per atom 
to be finite in the bulk limit, the exponent must exceed 
threshold values which depend on dimensionality of the 
lattice. For single dimensional systems, one must have 
a > = 1 and a > a| D = 2 for two dimensional 
crystals. 

We report on a calculation of the atomic root mean 
square deviation <5rms about positions of equilibrium in 
ID and 2D crystals. In section I, we discuss theoreti- 
cal methods used to calculate the partition function and 
hence calculate salient thermodynamic quantities by de- 
coupling the vibrational states used to gauge the integrity 
of long-range crystalline order such as 5rms- In Section 
II, we examine one dimensional systems, finding posi- 
tional order to be destroyed except in a long-range cou- 
pling scheme, (i.e. a power law dependence where the 
decay exponent a must lie between a^ D = 1 and an 
upper bound exponent a]P — 1.615 ± 0.005). In sec- 
tion III, we perform a similar analysis for two dimen- 
sional square lattices, where we generalize to an extended 
scheme, and find a gradual destruction of crystalline or- 
der with increasing system size L for short-ranged cou- 
plings. However, we find that a power law decay profile 
where a < = 3.15 ± 0.025 is sufficient to maintain 
positional order at finite temperatures. In addition to 
the square lattice, in Section IV we also examine tri- 
angular and honeycomb lattices, finding the ability of a 



long-range interaction between atoms to preserve crys- 
talline order is not affected by the specific type of lattice 
geometry under consideration, and the threshold expo- 
nents in all three cases are identical within the bounds 
of error in the calculations. Finally, in section V, we 
consider motion transverse to the plane of the crystal 
lattice for locally stiff dual layer systems where even if 
interaction between particles is taken to be long-ranged, 
we find the perpendicular motions rapidly compromise 
long-range order as the size of the system is increased. 

A salient component of our treatment is the ex- 
plicit accounting for atomic motions in discrete sys- 
tems. The harmonic approximation, which neglects an- 
harmonic terms in the potential set up by geometric ef- 
fects has been tested directly in the context of Monte 
Carlo simulations and found to be accurate to within 
one part in 10 3 for the systems considered-. 

Since our interest is in equilibrium thermodynamic 
characteristics of the system, we begin with the lattice 
potential 

N m 

i=i j=i 

where rii is the number of neighbors corresponding to 
each atom in the system indexed with the label i, and N 
is the total number of particles contained in the crystal 
lattice. The factor of 1/2 in the lattice potential ex- 
pression is present to compensate for double counting of 
the energy associated with individual "bonds" between 
atomic pairs i and j. 

For small deviations from equilibrium positions, a 
"harmonic approximation" is possible, and one finds in- 
stead 

V=\Y.f.^-l%)\ (2) 

i=i j=i 

the first nonzero term of a Taylor expansion of V(rij) 
where Uj ~ and Kij is the second derivative of the 
potential V(rij) at r.^ = l^. In the results we report on 
here, we restrict attention to temperatures below those 
which would cause melting in the bulk (determined by 
the Lindemann criterion), where the harmonic approxi- 
mation would tend to fare well. A primary issue of in- 
terest is whether there is any finite temperature range 
where long-range positional order is intact, and our cal- 
culations are in the context of temperatures not of the 
magnitude that would disrupt the bonding topology and 
create dislocations, but thermal regimes considerably be- 
low the temperature range which might begin to rupture 
bonds between neighbors. For covalent solids such as two 
dimensional sheets of graphene and carbon nanotubes 
where energies stored in covalent bonds are far in excess 
of k^T at 300K, even room temperature may be consid- 
ered a "low" temperature in the sense of being consider- 
ably below temperature scales where thermal fluctuations 
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would perturb the local bonding scheme in a significant 
way. 

In one dimension, the bonds are collinear, and the po- 
tential will remain quadratic as the energies of all bonds 
between atoms are summed. However, in two dimen- 
sional geometries, restoring forces to oppose displace- 
ments from equilibrium will be exerted in different di- 
rections along distinct bond axes between an interacting 
pair of atoms. As a consequence, it will be necessary 
to make an additional harmonic approximation in order 
to obtain a quadratic Hamiltonian and subsequently ex- 
ploit translational invariance for the regular lattices we 
examine. 

In general, a bond length Zy between atoms i and j 
will appear as 
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where A?f 



A 0y 



and A°/ 



(zf — Zj). Thus, the potential energy stored in the bonds 
depends only on the difference of coordinates such as, 
e.g., for the equilibrium x coordinate differences and 
(Sf — SJ) for differences constructed from the correspond- 
ing shifts from equilibrium. If the latter are sufficiently 
small in relation to the former, it is appropriate to ex- 

and to quadratic order one 
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pand about A?? 



A ay , 



A 0z , 



where Ay is a 



will have (Zy - Z°) 2 ps Ay • (Si - Sj 

unit vector formed by subtracting the position vectors x? 
and x° corresponding to the atom i and its neighbor j, 

such that Ay = (x° — x°)/||x° — x°||. Hence, the atomic 
potential may be written to quadratic order in Si and Sj , 
and one has in particular 
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Ay • (Si - Sj) 



(4) 



It will be necessary to solve an eigenvalue problem to 
decouple the vibrational modes. However, with Fourier 
analysis, the problem may be reduced to the task of di- 
agonalizing a 2 x 2 matrix, a 4 x 4 matrix in the case of 
a lattice with a honeycomb geometry, or at most a 6 x 6 
matrix for the case of the locally stiff dual-layer system, 
even in cases where the coupling scheme is extended to 
encompass many neighbors for each atomic member. We 
will consider systems in one dimension where we show 
that long-range crystalline order may be stabilized in the 
case of slowly decaying power law potentials, but not for 
localized exponentially decaying coupling schemes. We 
also examine two dimensional lattice geometries, and find 
similar phenomena; again, a long-range power law decay 
is needed to preserve crystalline order at finite tempera- 
tures. 

Finally, we are careful to restrict motion to collinear 
displacements in the context of ID systems and intrapla- 
nar motion for the two dimensional crystals. The lattices 



are very easily disturbed by transverse displacements, 
and we find that relaxing the collinear and coplanar re- 
strictions in ID and 2D yields mean square fluctuations 
which diverge rapidly with increasing system size N. As 
we find with explicit calculation in section V, this rapid 
(i.e. at a linear) growth in N occurs even if the interac- 
tion between atoms is a long-ranged power law decay in 
locally stiff dual layer crystal geometries. 

We use the results for the eigenvalues for the vibra- 
tional modes to calculate thermodynamic properties re- 
lated to crystalline order such as the thermally aver- 
aged mean square fluctuations about equilibrium per site, 
$rms- As noted elsewhere^ and summarized here, the 
mean square displacements about equilibrium may be 
calculated in terms of the eigenvalues for the vibrational 
states. 

In terms of the vibrational modes, the lattice energy 
may be written as 



1 3M 



(5) 



with M the total number of particles contained in the 
lattice. The connection between the vibrational states 
and the mean square fluctuations is 

1 M 

(<W 2 = £<(^) 2 + W + (W 2 ) (6) 

i=l 

With the eigenvectors indexed with the label a and using, 
e.g., Sf = X)a*fi c aV'a, to express the displacements in 
terms of the eigenvectors, one finds for a specific system 
configuration 
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s 2 = 



i—l a — 1 a ' —i 



(7) 



In the thermal average, the factor c a c a > will be as often 
positive as negative, and hence only in the case a — a 
will there be a net contribution to the thermal average 
(<5rms) 2 - If the vectors are taken to be normalized, one 
finds 



RMS/ 



1 3M 



(S) 



Q = l 



With the lattice energy expressed in this way, the parti- 
tion function becomes of a product of Gaussian integrals, 
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dc a = J , 
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1/2 



(9) 



where r = k^T. Finally, a thermal derivative of Z leads 
to 



3 A/ 



(W) 2 -^Ln(Z) = ^A- 1 r 



(10) 
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Hence, for the mean square deviation, we have 



\ 



3M 



RMS 



(11) 



= 1 



The term in the radical is not temperature dependent, 
but is instead determined by characteristics of the lattice 
geometry and the bonding scheme between atomic mem- 
bers. In this work, we calculate <5^ M g, a mean square 
RMS deviation normalized with respect to temperature. 
Zero eigenvalues are artifacts of the periodic boundary 
conditions, correspond to global translations of the lat- 
tice, and are excluded from the sum. 
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FIG. 1: Mean square deviation curves are graphed versus 
system size N in inset (a). The vertical axis of the graph in 
inset (b) is (5rms) 2 to emphasize the linear behavior of the 
square of the mean square fluctuations for a broader range of 
system sizes. 



II. SYSTEMS IN ONE DIMENSION 

In the ID systems we consider, only longitudinal dis- 
placements are examined; similarly, for the two dimen- 
sional geometries, lattice vibrations are considered to be 
confined to the two dimensional plane with no transverse 
motions considered. Periodic boundary conditions are 
implemented in both the one and two dimensional cases. 
An important characteristic of systems in one dimension 
is the fact that all bonds are collinear, and hence there 
is no purely geometric source of anharmonic effects. The 
lattice potential energy will have the form 



(12) 



1=1 m=l 



Since only longitudinal motions are considered, the label 
"x" that would normally appear as a subscript on the "<5" 
symbols is suppressed. The sum recorded in Eq.[12]is con- 
figured to avoid the redundant summation over bonds, 
and the counting "1/2" factor is not required. To maxi- 
mize the number of neighbors coupled to any particular 
atom while avoiding multiple couplings to the same atom 
via the periodicity condition, we set n — (A ~ l)/2 and 
we always consider an odd number of atomic members. 

It is convenient to operate in terms of Fourier compo- 
nents, where we have Si = e l/c '4; on substitution, the 
expression for the lattice energy has the form 



2J K m (l - cos mk x ) 
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(13) 



which has been diagonalized with the use of Fourier 
components Sk- The eigenvalues are given by Afc x = 



2 y ^K m (l — cos km), and the normalized thermally in- 
has the form 

1/2 



m— 1 

duced shift S 



RMS 



J RMS 




(14) 



We first examine a localized potential which in the ID 
context would certainly be expected to yield a divergent 
mean square fluctuation ^ MS with increasing system size 
A. As a companion result to gain complementary insight, 
we also calculate the density of states corresponding to 
the system. One merit of obtaining the vibrational den- 
sity of states is the fact that it may be computed in the 
thermodynamic limit without encountering divergences 
with the aid of Monte Carlo sampling. On the other 
hand, the divergence or convergence of the mean square 
deviations will be signaled by specific signatures in the 
low eigenvalue regime of the density of states without the 
need for an extrapolation to the thermodynamic limit. 

In the case where interactions are confined to near- 
est neighbors, the eigenvalues Afc^. have the simple form 
2Kq{1 — cos k x x). The corresponding thermally averaged 
^rms va hres m the case of finite systems are shown in 
the graph in Fig. [TJ and there is a steady rise with in 
the RMS fluctuations with increasing n. The expansion 
of the mean square deviations from equilibrium is sub- 
linear, but it may be shown that the increase continues 
indefinitely (i.e. diverges in the thermodynamic limit) 
by graphing instead (£>rms) 2 ' as m the inset of Fig. Q] 
The dependence on A quickly reverts to an asymptoti- 
cally linear increase with A, and to a good approximation 
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RMS 



A 1 / 2 for moderate to large systems. 



We next extend the coupling scheme to many neigh- 
bors where the coupling decays at an exponential rate, as 
might be found at least on a qualitative level for a cova- 
lently bonded system where the rapidly decaying overlap 
of the orbitals of atomic neighbors (and hence the mag- 
nitude of the exchange coupling) has an asymptotically 
exponential decay as the separation between the pair of 
atoms becomes sufficiently large. The lattice energy will 
have the form 



JV 



'H+mj 



(15) 



1=1 m=l 



Hence in terms of Fourier components, the total lattice 
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potential becomes 



v c i P = J2 2K ° 



[1 — cosmfc] 



(16) 



Again, operating in terms of Fourier components decou- 
ples the modes, and the appropriate eigenvalues A& are 
given by 



A* 



-7m 



[1 — cos mk] 



(17) 



where the prefactor K has been suppressed. In addi- 
tion, the label u x" on the wave vector has also been sup- 
pressed for the sake of convenience. Although the cou- 
pling scheme is extended to many neighbors, the poten- 
tial is in an important sense still a local interaction due 
to its rapid decay, where the appropriate length scale is 
the inverse decay rate 7 -1 . By appealing to the formula 
for a geometric sum, r + r 2 + ... + r n — (r — r n+1 )/(l — r) 
(where r is taken to be complex and \r\ < 1) and using 



the fact that cos mk = i(e* 



c ), one may obtain 



an explicit expression for A& which does not require the 
intermediate summation. Applying the geometric series 
formula for a finite series yields 



At 
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(18) 
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Combining the last two fractional terms gives 

( cos k + e _7< - n+1 ' cos nk \ 
-e~ 7 — e~ 7 ™ cos(n + l)k] 



e 7 - 1 



V 



cosh 7 — cos k 



J 



(19) 



a tidier and computationally convenient expression to use 
in calculating the RMS displacements <5r.ms- 

As n becomes large, terms proportional to e -7 ™ 
quickly become suppressed by the rapid exponential de- 
cay. Hence, for n 7 -1 , one will obtain 



A* 



cos(fe) 



e 7 — 1 cosh(7) — cos(fc) 



(20) 



The results are shown in Fig. [3] for the scaling of <52 MS 
with respect to the size N of the system. To keep the re- 
sults for different values of 7 on the same footing, we use 
the prefactor Kq as a normalization of the coupling with 
Kq = e 7 . A similar procedure is also used in calculations 
involving exponentially decaying extended couplings in 
2D. In this manner, the convergence to the results for 
the case where only nearest neighbors are involved in the 
coupling scheme is easier to see. 

In the main part of the graph, the square of the RMS 
deviation is graphed with respect to system size N for a 



broad range of system sizes. The curves corresponding 
to the different decay constants are asymptotically linear 
in the system size although the slopes decrease with de- 
creasing 7 as the coupling becomes longer in range. The 
inset of the graph shows a closer view of ^rms- Each of 
the curves rises steadily for sufficiently large N, notwith- 
standing non-monotonicities for small to moderate N in 
the case 7 = 0.25 where the decay of the interaction is rel- 
atively slow. In latter case, there is competition between 
thermal fluctuations and an increase in lattice rigidity 
which occurs as the linear crystal grows, providing atoms 
with more neighbors. Eventually, however, N exceeds the 
length scale 7 -1 of the coupling between atoms, and the 
balance shifts in favor of thermal fluctuations. The lat- 
ter increase in importance with increasing N and thus 
eventually destroy long-range crystalline order. 

We also examine the eigenvalue density of states for 
different decay rates 7. With the range of the potential 
being set by 7 -1 , larger values of 7 would correspond to 
a more rapid decay and a shorter range of the interaction 
between neighbors. In calculating the density of states, 
we use a Monte Carlo sampling process where the val- 
ues of k are not quantized, permitting one to genuinely 
achieve the bulk limit for the purpose of obtaining the 
density of states. To obtain a smooth curve a large num- 
ber of eigenvalues (i.e. 2.5 x 10 8 for the histograms corre- 
sponding to the exponentially decaying coupling scheme) 
are sampled. The formula given in the continuum limit 
in Eq. [5D] is the appropriate expression to use for Afc in 
the Monte Carlo sampling process. 

The normalized density of states for a range of decay 
constants 7 appears in Fig. [3] Panel (a) is a standard 
plot with the density of states on the vertical axis, while 
to facilitate the viewing of the DOS curves, the logarithm 
of the DOS curves in shown in panel (b). Even for rel- 
atively long-range cases such as 7 = 0.25, the density of 
states retains the "U" -shaped profile of the nearest neigh- 
bor case. The latter corresponds effectively to 7 = 00, 
and is shown in red in the graphs. For convenience in 
comparing results, we again choose Kq = e 7 in calculat- 
ing the DOS curves. The convergence to the 7 = 00 case 
with increasing 7 is evident for the case 7 = 2.0, where 
close agreement with the DOS calculated for the nearest 
neighbor case is evident in panel (b) of Fig. [T] In the 
latter, the ordinate is chosen to be log 10 (DOS) to help 
show the structure of the density of states curves more 
clearly. 

Finally, we examine a genuinely long-range inter- 
atomic potential with a power law decay profile. The 
lattice potential energy will have the form 



(21) 



1=1 m=l 



where a is the decay exponent of the power law inter- 
action (a > 1), and again n — (N — l)/2. In terms of 
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FIG. 2: (Color Online) Mean square deviation curves are 
shown for a range of 7 values. The vertical axis of the 
main graph is (<5r,ms) 2 to emphasize the linear behavior of 
the square of the mean square fluctuations, while the inset 
is a graph of <5rms with respect to N for a smaller range of 
system sizes. 2.5 x 10 s eigenvalues were sampled to calculate 
the DOS curves. 
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FIG. 3: (Color Online) Eigenvalue Density of States (DOS) 
for 7 values ranging from 7 = 0.25 to 7 = 2.0 to, effectively, 
7 = 00 with the DOS on the ordinate for panel (b). Panel 
(b) shows the density of states as well, but the ordinate is 
log 10 (DOS). 



Fourier components, one will have 



E = J2Y1 2K om~ a (l - cosmk)\5 k \ 



(22) 



fe m=l 



Hence, the modes are now decoupled with eigenvalues 

n 

given by Afc = m~ a (l — cosmfc); we evaluate this ex- 

m— 1 

pression directly in order to obtain <5r M s and the DOS 
profile appropriate to particular exponent a in the bulk 
limit. Again, we calculate S^ MS and generate plots with 
respect to system size. To test for divergence or con- 
vergence in the bulk limit, it is useful also to prepare 
log- log plots (we use base ten logarithms in all cases), 
and the results appear in the inset of Figure 21 We ex- 
amine systems ranging in size from N = 3 to on the order 
of a few hundred thousand atomic members. A crucial 



question is whether there is a threshold value a c above 
a = 1 (where the lattice energy may diverge with in- 
creasing system size) below which long range crystalline 
order is stable with respect to thermal fluctuations in one 
dimensional lattices. 

To identify a]P , we calculate the normalized mean 
square fluctuations <5r M s with respect to system size n, 
producing log-log graphs. The highest value of a where 
the mean square deviations converge is identified as a]P , 
the upper limit for the decay exponent in the extended 
power law decay scheme where long-range crystalline or- 
der is still supported at finite temperatures. 

The mean square deviations, useful thermodynamic 
quantities with which to diagnose the presence or absence 
of long-range crystalline order, are shown in Fig. [4] and 
Fig. [5] (with the abscissa shown as a base ten logarithm 
over five decades of system sizes N). In Fig. [U 5^ MS 
curves are shown for a relatively wide range of a values. 
Over the broad range of systems on the horizontal axis, 
five orders of magnitude, the mean square displacements 
rise monotonically for a = 2.0 and a = 1.75, while <5r MS 
decreases steadily for a — 1.5 and a = 1.25. The curves 
suggest a decay exponent a c in the vicinity of a — 1.6 as a 
boundary between crystals where long-range order is un- 
stable at finite temperatures, and one dimensional solid 
where crystalline order is retained for T > 0. The inset 
is the corresponding log-log graph of the mean square 
fluctuations plotted for the same a values as in the main 
graph, which is a semi- logarithmic plot. 

In Fig. RMS deviation curves are shown for a 
tighter span of power law decay exponents (ranging from 
a = 1.60 to a = 1.65) to identify with greater accuracy 
the numerical value of a c . To facilitate the determina- 
tion of the exponent separating crystals with long-range 
order and those disrupted by thermal fluctuations, we 
place dark circles over the maxima of the S^ MS curves. 
For a > 1.625, the maxima are located at the edge of 
the plot, consistent with a steady increase (and likely di- 
vergence in the bulk limit) of the RMS curves. On the 
other hand, for a < 1.6125, the thermally averaged RMS 
deviations are non-monotonic, reaching a maximum for 
finite values of N and then declining, presumably to- 
ward a stable bulk value. We identify the boundary as 
a c = 1.615 ± 0.005. It should be emphasized that while 
long-range order is not supported for decay exponents in 
excess of cv c , the divergence of 5^ MS with increasing sys- 
tem size is nonetheless quite slow, sublinear in log 10 (AT), 
whereas a strictly logarithmically diverging mean square 
deviation would instead rise at a more rapid linear rate. 

To obtain information complementary to the RMS 
fluctuations, we again calculate the eigenvalue density 
of states. We also use Monte Carlo sampling where wave 
numbers are chosen at random, with uniform probability, 
to calculate the vibrational density of states with the re- 
sults shown in Fig. [6] The double sum in Eq. [22] requires 
careful consideration, in that one must be certain that 
enough terms have been included in the inner sum that a 
convergent result is obtained. To be certain convergence 
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FIG. 4: (Color Online) <5r,ms i s plotted on a semi-logarithmic 
(all logs are base ten) scale in the main draft for selected 
values of a for one dimensional systems. The inset is a log- 
log graph of the mean square deviations for the same values 
of a as in the main graph. 
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FIG. 5: (Color Online) <5r,ms i s plotted on a semi- logarithmic 
scale for a relatively tight range of a values. The large dark 
circles indicate maxima in the mean square deviation curves; 
it is concluded that a c = 1.615 ± 0.005. 



has been achieved, we prepare eigenvalue histograms for 
successive doublings of the number of terms contained in 
the inner sum indexed by m. The number of terms which 
must be included in order to attain suitable convergence 
increases with decreasing a for crystal lattices where the 
coupling is more slowly decaying. In general, however, 
the oscillatory cosine term in Eq. does act to somewhat 
hasten convergence and hence limit the number of terms 
which need to be summed. 
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FIG. 6: (Color Online) Eigenvalue Density of States (DOS) 
for the one dimensional lattice with a power law decay where 
the decay exponent a varies from a = 1.5 in panel (a) to 
a = 2.5 in panel (d). 
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FIG. 7: (Color Online) Three distinct lattice geometries are 
shown. Panel (a) is a portion of a square lattice, panel (b) rep- 
resents the triangular lattice, and panel (c) shows the honey- 
comb lattice with inequivalent sites represented with distinct 
colors. 



III. TWO DIMENSIONAL CRYSTALS 

For the case of a two dimensional system, the analysis 
is in many respects parallel to that applied for the one 
dimensional lattices. However, the additional dimension 
makes available richer choices for the lattice geometry. 
We examine various coupling schemes for three types of 
lattices; the square lattice, the triangular lattice, and the 
honeycomb lattice. In Fig. [7] panel (a) represents the 
square lattice, the triangular lattice is depicted in panel 
(b), and the honeycomb lattice appears in panel (c). A 
peculiarity of the honeycomb lattice is the presence of in- 
equivalent sites, and this characteristic is highlighted in 
panel (c) of Fig. [7] where different colors are used in label- 
ing the sites. Although the geometries we examine have 
different characteristics, the essential qualitative charac- 
teristics and the most salient physics are found to share 
much in common. 

We first consider the square lattice, and we ini- 
tially take into account only interactions between nearest 
neighbors where at the present level of approximation the 
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lattice lacks rigidity. The lattice energy is 
N-l 



E 



K 

~2 ^ 



^[(^i-^) 2 + (^ +1 

i,j=0 



S v ) 2 ~\ 

ij+1 °ij> 1 



(23) 



We express the displacements in terms of Fourier compo- 
nents with, e.g., Sfj = '^^S'£e I ( k, * t+kv: >\ with / being the 

k 

imaginary unit. In terms of 5£. and <5^, the energy has 
the form 

2 ^ 

k 



[(1 - cosfc x )l<5kT + (1 - cosk y )\Sm (24) 



in this manner the degrees of freedom are decoupled. In- 
spection of Eq. [24] reveals that the eigenvalues are 2N 
fold degenerate and identical to the eigenvalues obtained 
for the case of the one dimensional crystal where only 
interactions between nearest neighbors were considered. 
Since the eigenvalues are the same as those in the ID case 
with interactions only between nearest neighbors, crys- 
talline order is readily disrupted by thermal fluctuations. 
Hence, (^yyis) 2 wm sc& l e with N just as was the case for 
the ID counterpart. 

If one takes into account coupling to next-nearest 
neighbors as well, then the lattice energy in real space 
is 



N-l 
i,j=0 



[x ■ (8 i+ ij 



ij+1 



V2 



H+lj-1 



2 \ 



(25) 



Operating in terms of Fourier components, one diagonal- 
izes the matrix 



2 — cos k x 
cos fe T cos k„ 



sin k T sin fc„ 



sin k x sin k y 

2 — cos k y 
- cos k T cos k„ 



The eigenvalues are given by 

A± = (4 — cos k x cos k y — 2 cos k x cos k y 
± J (cos k x — cos ky) 2 +4 sin 2 k x sin 2 k y 



(26) 



(27) 
(28) 



The results for the mean square deviations 6^ MS ap- 
pear in Fig. HI One may also calculate the vibrational 
DOS, and the results appear in panel (a) of Fig. [5J The 
introduction of next-nearest neighbor interactions is very 
effective in reducing the deleterious effect of thermal fluc- 
tuations on long-range crystalline order, though there is 
still a weak divergence in S^ MS in the bulk limit. The 
square (^ MS ) 2 quickly assumes an asymptotically linear 
form with respect to log 10 N. The much slower increase 
of the RMS deviations with N is reflected in the DOS pro- 
file, where instead of exhibiting a sharp cusp in the low 
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FIG. 8: (Color Online) The mean square displacements for 
the square lattice with nearest and next-nearest neighbor cou- 
plings. The main graph shows (5r,ms) 2 j with the correspond- 
ing vibrational density of states graphed in inset (a) and the 
raw RMS displacements <5rms displayed in inset (b). 



eigenvalue regime, the DOS curve terminates smoothly. 
However, the fact that the DOS tends to a finite value 
as the eigenvalue vanishes is still enough to cause a di- 
vergence in the mean square displacements from equilib- 
rium. 

We next examine a general case where there are inter- 
actions with many neighbors. In real space, the energy 
stored in the lattice has the form 



N-l n 

E E 

i,j=—n l,m=— 



\K{r lm ) {k lm ■ [5 i+lj+m ~ 4]) 2 (29) 



where the inner 



tin 



factor compensates for multiple 



counting of bond energies and the choice n = (N — l)/2 
allows each atomic member to interact with all of the 
atoms contained in crystal while avoiding multiple in- 
teractions with the same particle. Since A; TO = [Z,m], 
the appropriate unit vector directed between particles 

given the labels "i + l,j + m" and "if is A; TO = 
[/, ml/y/l 2 + m 2 . Again, we may decouple the vibrational 
modes by expressing the coordinate shifts in terms of 
Fourier components. The lattice potential energy may 
then be written as 



N-i n 

* = E ^ 

i,j=0 l.m= — 



(30) 



with the "1/2" factor present to compensate for redun- 
dant bond counting. The range radius r\ m is defined 
with r; m = y/l 2 + m 2 , with the full vector given by 
ri m = Ix + my. In terms of Fourier components, one 
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will have 



1 



£ = 5 E £ 



A-(r, m ) 



[1 - cos(k ■ r lm )] 



(31) 



' in? 



l 2 \S^ 2 



m 2 \S y ^ 2 



Hence, in order to to decouple the vibrational modes, one 
must diagonalize the matrix 



2 £ 



K(r lr 



[1 - cos(k • r ;m )] 



I 2 



mZ 

™2 



> (32) 



which may also be written as 



K(r lm )[l - cos(k-r im )] 

I m— — n 



Af Af Af Af 

/m Ira Ira ira 

Af Af Af Af 



, (33) 



a representation which will prove more compact for more 
complicated systems such as the honeycomb lattice crys- 
tals with more than one layer in the direction transverse 
to the crystal plane, examined in Section V. 

We first turn to the case of an exponentially decay- 
ing coupling scheme, and we calculate the S^ MS curves 
with respect to the size N of the system. The results for 
the thermally averaged means square displacements are 
shown in Fig. |H] for a range of different decay constants 
7. The computational burden of calculating the auxiliary 
sum will grow with N, but one aspect of the exponential 
decay that is of assistance in the calculations is the fact 
that the sum may be safely truncated when the distance 
between interacting atoms becomes several times greater 
than the range of the short-ranged coupling (i.e. terms 
beyond beyond m and I pairs such that \/P + m 2 3> 
need not be included. In particular, we obtain results 
which are very well converged if we discard terms beyond 
20 decay lengths 7 -1 . Ultimately, thermally induced de- 
viations from equilibrium destroy long-range order, and 
the RMS deviations diverge slowly [(S^ MS ) 2 again scales 
linearly with log 10 (7V)], but the rate of divergence de- 
creases with decreasing 7. In particular, as the range 
7 _1 of the inter-atomic coupling is increased, the slope of 
the graph of (S^ MS ) 2 with respect to log 10 (JV) decreases, 
although the RMS deviations eventually still diverge in 
the thermodynamic limit. 

The DOS is also calculated, with results appearing in 
Fig. [10] for a range of 7 values. We use Monte Carlo sam- 
pling to choose k x and k y from a continuum range, and 
thereby operate in the thermodynamic limit for the pur- 
pose of calculating the DOS curves. At least 10 6 eigenval- 
ues are sampled in generating the DOS curves. The low 
eigenvalue region of the DOS graph is very similar to the 
corresponding regime of the density of states where only 
interactions with nearest and next-nearest neighbors are 
included. 

Next, we examine the much slower power law decays 
Kij = Kr~, a in the inter-atomic separation ru where the 
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FIG. 9: (Color Online) Mean square deviations for an ex- 
ponentially decaying interaction. The main graph shows 
(<5rms) 2 > while the raw mean square deviations are plotted 
in the inset of the Figure. 




2.0 4.0 6.0 8.0 10.0 

Eigenvalue Magnitude 



0.45 


(d) | y = 0.5; square lattice | 




0.35 






0.25 








0.15 






0.05 







2.0 6.0 10.0 14.0 

Eigenvalue Magnitude 



0.0 5.0 10.0 15.0 20.0 25.0 

Eigenvalue Magnitude 



FIG. 10: (Color Online) Vibrational density of states for an 
exponentially decaying coupling. The DOS curves are plotted 
for assorted values of the decay parameter 7. 



exponent a controls the rate of the decay, long-ranged in 
the sense that there is not a length scale to set the range 
of the coupling. Again, we first calculate the mean square 
displacements with respect to the size of the system, and 
then we examine the density of states for the eigenvalues. 
The (Tr MS results are shown in Fig. [TO 

Computational subtleties similar to those encountered 
for the case of the one dimensional solid in with a long- 
ranged coupling scheme must be carefully navigated since 
the higher dimensionality {d — 2) will cause the compu- 
tational burden to grow even more rapidly (nominally as 
L 4 if interactions with all neighbors are included) with 
system size. Again, we use Monte Carlo sampling to se- 



10 




Log 10 (AO 



FIG. 11: (Color Online) Mean square deviations for an inter- 
action decaying as a power law for decay exponents a near 
the threshold o^P . The inset shows a broader view, while the 
main graph is a closer view of the transition from converging 
to diverging 5r MS curves. 
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FIG. 12: (Color Online) The vibrational density of states for 
the two dimensional square lattice with a power law decay 
interaction. DOS curves for various values of the decay pa- 
rameter a are shown, with traces for different values of the 
truncation length Na on the same plot. 



lect wavevectors and accumulate eigenvalues to build up 
the vibrational DOS. The auxiliary sum over dummy in- 
dices I and m giving the eigenvalue would in principle 
contain an infinite number of terms (in the bulk limit, 
each atom in the crystal would have an infinite number 
of neighbors), but we truncate the sum at a finite range. 
The presence of sinusoidal terms in the sum, as in the 
corresponding ID case, provides an oscillatory element 
and will hasten the convergence of the sum, thereby re- 
ducing the computational burden. We check convergence 
with respect to the truncation range by calculating the 
DOS with successive doublings of the truncation length 
TVa until the DOS profile ceases to change with addi- 
tional doublings of the system size. The results for the 
vibrational DOS appear in Fig. [12] One notes that the 
convergence with respect to the truncation radius is least 
rapid for a — 2.5. However, the graphs are relatively 
well converged for a = a 2D and higher values of the de- 
cay exponent. When a is in the vicinity of a 2D , there 
is very little support in the low eigenvalue regime. On 
the other hand, with increasing a the interaction decays 
more rapidly, and ultimately the histogram amplitude in 
the zero eigenvalue limit rises to a finite value, contribut- 
ing to a divergence in the mean square fluctuations with 
increasing system size. 



IV. ALTERNATE TWO DIMENSIONAL 
GEOMETRIES 

The treatment in the case of the triangular and honey- 
comb lattices is very similar to the approach used in the 
case of the square lattice. Interestingly, the triangular 
lattice is rigid with only a nearest neighbor interaction 



in the context of the harmonic approximation, and the 
mere inclusion of nearest neighbors is enough to set up 
quasi-long range order where thermally induced fluctua- 
tions about equilibrium diverge very slowly [i.e. (^rms) 2 
increases as the logarithm of the system size just as in 
short-ranged extended interactions for the square lattice] . 
For the triangular lattice, the lattice energy in real space 
has the form 



N-i 

£ = f E 

i,j=0 



( 



H+lj 



-Si 



{¥ + #y) • (4+i - 1 



[\x - ^-y) ■ [Si+ij-i - 5, 



\ 



I 



(34) 



and the eigenvalues for the decoupled vibrational modes 
are obtained by diagonalizing the 2x2 matrix 



3 — 2 cos k 

~ 2 

(cos[fc 



■ /c-^J cos ky 



v3 
2 



} cos k y 



"2 COS ky 



cos[ky k x ] 



yielding 

[3 — cos 



cos k y — cos(k y — k x )\ 



, I cos 2 k x + cos 2 ky + cos 2 (k y — k x ) — cos k x cos k y 
1 — cos k y cos(k y — k x ) — cos k x cos(k y ~k x ) 



(35) 



(36) 



The results for the mean square fluctuations as well as 
the vibrational DOS are given in Fig. [T31 

We generalize the nearest-neighbor case to an extended 
scheme where each atomic member may interact with 
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FIG. 13: (Color Online) The mean square deviation for the 
triangular lattice with nearest neighbor couplings. The main 
graph is a semi- logarithmic plot of <5rms , while inset (a) shows 
the vibrational density of states. A graph of the raw <5rms 
appears in panel (b). 
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FIG. 14: (Color Online) Normalized Mean square curves and 
eigenvalue density of states for the triangular lattice with an 
exponential coupling scheme. The main plot shows (5rms) 2 
with respect to the base ten logarithm of the system size N, 
while inset (a) is the corresponding graph for (S^Mg. Inset 
(b) shows the density states curve for the triangular lattice 
for 2.5 x 10 8 eigenvalues sampled in the bulk limit via Monte 
Carlo 



many neighbors. In real space, the lattice potential en- 
ergy may be written as 



E = i ^ ^ ~K(ri m )Ai m -(S i+ i tj+m - Sij), (37) 



iV 



i.j— l,m=—n 



pairs in the triangular lattice geometry. After expressing 
the displacements in terms of Fourier components, one 
calculates the eigenvalues of the 2x2 matrix 



9h 

l,m=—n 



Af Af 

lm i 

Af A? 

lm i 



Af Af ' 

lm lm 

Af Af 

lm lm. 



(38) 



where gi m = [1 — cos(k x l + k y m)]K{ri m ). As in the case 
of the square lattice, we consider for the triangular lat- 
tices an exponentially decaying coupling scheme, and the 
results are shown in Fig. [14] for over four decades of sys- 
tem sizes. Again, the quantity (S^ MS ) 2 increases linearly 
as /3 log 10 N with the slope j3 decreasing with decreasing 
decay rate 7 (and hence increasing range of the coupling). 

We also prepare graphs of the vibrational density of 
states, shown in the four graphs in Fig. [15] for a range 
of values of the decay constant 7. The wide separation 
between the RMS curves corresponding to 7 = 2.0 and 
7 = 1.0, 7 = 0.75, and 7 = 0.5 is mirrored in the DOS 
curves where for the smaller decay rates the eigenvalue 
histogram curve intersects the ordinate with very low am- 
plitudes. On the other hand, for the more rapid decay 
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where the components of the unit vector Ai ri 
(I + m/2)/rj m and Af m = V3m/2r ;m , where n m = 
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FIG. 15: (Color Online) Density of states curves for 7 = 2.0, 
7 = 1.0, 7 = 0.75, and 7 = 0.50 in panels (a), (b), (c), and (d) 
respectively for the triangular lattice with an exponentially 
decaying coupling scheme. 



where 7 = 2.0, the amplitude in the regime of low eigen- 
values is much higher, and the DOS graph resembles that 
of the nearest neighbor case to a much greater degree 
than DOS profiles corresponding to lower decay rates of 
the exponential coupling. 

As for the square lattice geometry, we examine a long- 
ranged power law interaction in the context of triangular 
lattices. The mean square deviations from equilibrium 
are graphed in Fig. 1161 with the inset of the plot showing 
a closer view of the (^rms) 2 curves. A salient question is 
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FIG. 16: (Color Online) The square of the normalized RMS 
deviations for various values of the decay exponent a, with 
the inset showing a closer view of the (<5rms) 2 curves. 



if lattices with geometries which differ from those of the 
square lattice will exhibit long-range crystalline order for 
the same range of decay exponents a as in the context 
of the square lattice. We find up to the bounds of error 
a^P = 3.15 ± 0.025 calculated for the triangular lattice 
to be identical to the threshold exponent for the square 
lattice. 

We show the corresponding eigenvalue histograms for 
the power law decay for the decay exponents a = 5.0, 
a = 4.0, a = 3.5, and a = 3.125 in panels (a), (b), 
(c), and (d) respectively of Fig. [171 While the eigenvalue 
histograms plotted in the panel (a) and panel (b) corre- 
spond to decay exponents significantly higher than a^P , 
the DOS curve in panel (c), is plotted for a decay ex- 
ponent only slightly above the threshold value, and the 
histogram in panel (d) of Fig. corresponds to a value of a 
just below (though very nearly equal) to a 20 . Whereas 
the DOS curves in panel (a) and (b) clearly tend to a 
finite value as the eigenvalue approaches zero, the ampli- 
tude for the slower decay a = 3.125 and tends to zero 
in the the limit that the eigenvalue is very small; for the 
case a = 3.5, the amplitude reaches a finite but very 
small value in the zero eigenvalue limit. A DOS ampli- 
tude tending to zero in the low eigenvalue limit, as seen 
for a — 3.125 is consistent with the preservation of long- 
range crystalline order indicated in the convergence of 
the mean square deviations graphed in Fig. 1161 

We examine the honeycomb lattice, which differs from 
the square in triangular lattices in that it possesses two 
inequivalent sites (labeled "A" and "B" for convenience). 
We again appeal to translational invariance, operating in 
terms of Fourier components, to decouple the vibrational 
modes for the honeycomb lattice. The relationship of 
sites of type "A" and "B" to nearest neighbors is illus- 
trated in Fig. [201 

Following this labeling convention, the lattice energy 
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FIG. 17: (Color Online) Eigenvalue histogram curves plotted 
for a = 5.0, a — 4.0, a = 3.5, and a — 3.125 respectively. The 
red trace corresponds to a relatively short truncation radius 
where Na = 25, and the truncation is less drastic for the blue 
DOS profiles where N A = 50. 
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FIG. 18: (Color Online) Labeling and indexing scheme for 
the honeycomb lattice for inequivalent sites labeled "A" and 
"B" , and their immediate vicinity. 



in real space has the form 



N-l 
i,3=0 



V 2 

< I ( — 



2 ^ 



s A - s 3 



(39) 



where it is sufficient to sum over the three bonds sur- 
rounding the atoms labeled "A" with no factor of | 
needed to compensate for double counting. In Fourier 
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space, the energy stored in the lattice has the form 



E = 



K 



E 



(l<ri 2 + I^T + iei 2 + lCf) ' 



■id 



_i(l + e ik *)e- ik v(3S* x *5£ v + 5® y *5£ x 



4 



(40) 



In addition to Fourier decomposition, the diagonalization 
of a 4 x 4 matrix will be necessary to completely decou- 
ple the vibrational modes appropriate to the honeycomb 
lattice with the nearest neighbor coupling scheme; the 
matrix in question is 



c AxAx c AxAy c AxBx c AxBy 

c AyAx CAyAy c AyBx c AyBy 

CBxAx CBxAy CBxBx CBxBy 

CByAx CByAy CByBx CByBy 



A B 
St A 



(41) 



where A and B, and W is the Hermitian conjugate of B. 
The sub-matrices A and B are given by 



A = 



and 



B 



-1(1 + e" 



ik x \ 
-1) 



V3 ^g — ifex 



1) 



■y _ 1(1 



(42) 



(43) 



However, to obtain a crystal which is locally stiff, one 
must examine an extended scheme where each atomic 
member interactions with many neighbors. In real space, 
the lattice energy may be expressed as 



N 

E =J2 

i,j=0 



(44) 



( \K{r\ m ) 



E 



(A 



(A 



(t)x A 

lm x 



-A^y).(^ 



lm 



i-\-lj-\-m ij ■ 



i+lj+m ij I 



K(rt) (A 



lm 



+lj+m w ij 



J 



where the first and second terms in the sum take into ac- 
count interactions between atoms labeled "A" and "B", 
respectively; the identical form of "A- A" and "B-B" in- 
teraction terms is due to the fact that the "A" and "B" 
species both define triangular lattices, as illustrated in 
Fig. We take the lattice constant to be unity, and the 

components of the unit vector A\ m appropriate to the 

triangular sublattices are A^ = V3(l + m/2)/rj m and 

(3[Z 2 +m 2 + lm\) 1 ' 2 . 



Af« = (3m/2)/r* 



* m , where r*„ 




FIG. 19: (Color Online) Honeycomb lattice geometry showing 
the interpenetrating triangular lattices defined by the inequiv- 
alent sites labeled "A" and "B" . 



On the other hand, the components of Af^ used in calcu- 
lating interactions between "A" and "B" atoms are given 



by A; 

where 



:(ab) 

„ab 
lm 



-- (3Z 2 - 



■m/2)/r ab and A 
- 3m 2 + 3/? 



y(ab) 
lm 

3m - 



1) 



1/2 



ab 
lm' 



The exploitation of translational invariance by express- 
ing the displacements in terms of Fourier components 
reduces the decoupling of the vibrational modes to the 

A B 

diagonalization of the 4x4 matrix ' 



sub-matrices are given by 



^= E 



mi, 



l,m= — n 



lm lm 

Af (t) Af (t) 



St A 



lm lm 
lm lm 



, where the 



(45) 
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for A and 
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y(ab) Xx(ab) Xy(ab) Xy(ab) 



x lr 



(46) 



where g tm = (1 - cos[k x l + k v m])K(rj m ), and g 00 = 
to exclude self-interactions in the "A- A" and the "B-B" 
coupled pairs. 

As we did for the square and triangular lattices, we 
examine a short-ranged exponential interaction between 
atoms in the lattice with a length scale given by 7 _1 . On 
the other hand, as for the preceding two lattice geome- 
tries, we also consider a long-ranged power law decay 
with K(r) — Kor~ a . The <5r MS results for the mean 
square deviations for the exponential decay scheme are 
shown in Fig. [501 for a range of 7 values. As in the 
cases of the square and triangular lattices in the extended 
schemes, increasing the range 7 -1 of the coupling slows 
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Log 10 W 

FIG. 20: (Color Online) Mean square deviations for the hon- 
eycomb lattice for an extended coupling scheme where the 
interaction decays exponentially. The main graph is a semi- 
logarithmic plot of (<5rms) 2 f° r various values of the decay 
constant 7, while the inset is the corresponding plot of 8™ m . 



(but does not halt) the rate of divergence of the mean 
square fluctuations. The large separation between the 
RMS curves corresponding to 7 = 2.0 and the slower de- 
cays 7 = 1.0, 7 = 0.75, and 7 = 0.50 is consistent with 
changes in the density of states curves where the his- 
togram amplitude in the low eigenvalue regime is sharply 
diminished as 7 decreases from 7 = 2.0 to 7 = 1.0. 

In the case of a power law decay, results for the mean 
square deviations arc shown in the semi-logarithmic 
graphs in Fig. l23l where the systems sizes considered span 
two decades, with a closer view for a more restricted 
range of the decay exponent a in the main graph; the 
inset is a graph of RMS curves for a broader set of a val- 
ues. As in the cases of the square and triangular lattices, 
the threshold exponent a™ is determined by examin- 
ing whether the RMS curves converge or diverge for very 
large system sizes. In agreement with the square and tri- 
angular geometries, we find oqP = 3.15 ± 0.025 for the 
critical decay exponent. 

The eigenvalue histogram curves shown in Fig. 1231 are 
consistent with the behavior of the mean square deviation 
curves given in Fig. 1221 For the more rapid decays a — 
4.0 and a = 4.5, the amplitude of the density of states 
is finite, which eventually contributes to a divergence in 
^rms ■ The divergence in the mean square fluctuations is 
much slower for a = 3.5, a characteristic which is echoed 
in the eigenvalue histogram in panel (c) of Fig. 1231 where 
in the zero eigenvalue limit the histogram amplitude is 
finite but very small. Finally, for a = 3.125, just below 
a^ D , the density of states curve tends to zero in the low 
eigenvalue limit. 
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FIG. 21: (Color Online) Density of states curves in the case 
of the honeycomb lattice for 7 = 2.0, 7 = 1.0, 7 = 0.75, and 
7 = 0.50 in panels (a), (b), (c), and (d) respectively. 




Log 10 (A0 

FIG. 22: (Color Online) Mean square deviations for long- 
range interactions with a power law decay exponent a for 
the honeycomb lattice geometry. The inset shows a relatively 
broad range of a values, whereas the main graph is a closer 
view. 



V. TRANSVERSE DISPLACEMENTS IN AN 
EXTENDED SCHEME 

In each of the preceding cases discussed in this work, 
thermally induced deviations of the lattice sites have 
been confined either to motion withing the lattice plane 
in the context of two dimensional systems, or collinear 
movements for the lattice in ID. However, we also con- 
sider displacements perpendicular to the lattice for two 
dimensional systems. To provide local stiffness of co- 
valent two dimensional lattices realized in nature where 
the finite thickness would provide rigidity with respect to 
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FIG. 23: (Color Online) Density of states curves in the case 
of the honeycomb lattice with a power law coupling scheme 
for a = 4.5, a = 4.0, a = 3.5, and a = 3.125 in panels (a), 
(b), (c), and (d) respectively. 



transverse displacements tending to push atoms above or 
below the plane of the crystal, we examine a dual layer 
geometry where an extended scheme provides a local stiff- 
ness. 

In fact, to provide as much stability as possible, we 
consider a long-range coupling scheme in the interaction 
between atoms within a layer as well as between layers 
decreases as a power law (with the decay exponent des- 
ignated a) in the separation between atomic species. 

The potential energy stored in the lattice is similar 
in abstract form to the corresponding expression for the 
honeycomb lattices, and is given by 



E = 
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E E 

l,7n— — n 
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A ab . (X B ^ xa 

^Im \ u i+lj+m 



\ 



S A ) 



(47) 



in a real space representation, where the unit vectors 

corresponding the intraplanar couplings A^ m have x and 
y components given by 



of Fourier components reduces the decoupling of the vi- 
brational modes to the diagonalization of a 6 x 6 matrix 
of the form where A and B are 3x3 sub-matrices, with 
W the Hermitian conjugate of B. The sub-matrices have 
the form 



^=E 



l.m——n 



/ jxx(ab) 
\ a im 
f ij/x(ab) 
\ a lm ™lm 
jzx(ah 

d i m 



,xx(s)^ 
L lm J 
M x { s )\ 



, ixy(ab) 
\ U lm 



1 lr, 



d 



'W(ab) 

lm ~ r tl ?m 
,zj/(ab 
a lm 



yy{s)\ 



a im 
,2/2 (ab) 

a lm 



lm 



(49) 



for A, where for the sake of brevity we have used the no- 



,xy(ah) _ 



tation, e.g. u, im 
matrix B is given by 



and d 
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lm 



The complex sub- 
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(50) 



Results for the mean square deviation <5p MS are shown 
in Fig. [M]for a variety of decay exponents ranging from 
a weak decay a = —2.5 to a considerably more rapid 
decay with the separation distance a = —6.0. A salient 
feature in each of the RMS curves shown in the graph 
is an asymptotic linear growth with the size N of the 
system, though the slope of the linear diverge in sys- 
tem size decreases with decreasing a. Broadly speaking, 
there are two regimes for each value of a in the varia- 
tion of the mean square deviation with system size. For 
relatively small system sizes, <5rms changes quite slowly 
with increasing N. However, ultimately the RMS fluc- 
tuations begin to grow more rapidly and eventually The 
size of the plateau where the mean square fluctuations ex- 
pand slowly is slightly broader for small values of the de- 
cay exponent a, and somewhat abbreviated for the more 
rapidly decaying coupling where a = —6. The latter 
is closer to what one would find in covalently bonded 
(but non-polar) systems such as graphene where Van der 
Waals interactions decreasing with the sixth power of 
the inter- atomic separation constitute the main source 
of long-range attraction between particles. Hence, Lon- 
don interactions would not prevent atomic displacements 
transverse to the lattice plane from destroying crystalline 
order. 



VI. CONCLUSIONS 



^Im — 



I 



Af (s) = 



(48) 



where the separation between interacting sites within a 
plane is rf m = Vl 2 + m 2 The components of the unit 

vector are identical to those of the planar case with 
the exception of a nonzero z component A^ b ' = 1 /rf r 
where the distance between sites in the two distinct lat- 
tice planes is rfj? = \Jvn? + I 2 + 1. Operating in terms 



,.ab 



We have examined the effect of thermally induced lat- 
tice vibrations on long range order in one and two dimen- 
sional crystals. In the case of crystals in one dimension, 
long-range positional order is (as expected) disrupted by 
thermal fluctuations with the 5^MS sca l m g as N 1 / 2 . For 
inherently long-ranged interactions scaling as power laws 
in the distance between interacting atoms, the divergence 
is much slower for a > — 1.615, while crystalline or- 
der is intact even at finite temperatures if a < a^P . 
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FIG. 24: (Color Online) Mean square deviations for long- 
range interactions with a power law decay exponent a. 

For two dimensional crystals, we find the same essen- 
tial phenomena with respect to thermodynamic stability 
of the crystal for square, triangular, and honeycomb lat- 
tices. For the latter two, thermal fluctuations destroy 
long-range crystalline order at finite temperatures, but 
the divergence in <5^ MS occurs very slowly with increas- 
ing system size. On the other hand, RMS deviations de- 
cay rapidly for simple square lattices where only nearest 
neighbor couplings are active. However, an extended cou- 
pling scheme to both nearest and next-nearest neighbors 
considerably mitigates the effect of thermal fluctuations 
on crystalline order, and the resulting slow divergence 



is quantitatively similar to that seen in the triangular 
and honeycomb lattices where the coupling is confined 
entirely to nearest-neighbors. 

When we extend the coupling to many neighbors, but 
still implement a short-ranged coupling (e.g. an exponen- 
tial decay with the range set by the inverse decay rate 
7 _1 ), we find qualitatively the same results as for the 
square lattice with couplings both to nearest and next- 
nearest neighbors, as well as the triangular and hexago- 
nal lattices in atomic members only interact with nearest 
neighbors with (<5rms) 2 scaling linearly with log 10 (iV). 
However, the slope of the linear dependence becomes 
smaller as 7 is decreased, as a longer-ranged coupling is 
more effective at suppressing the effects of thermal fluc- 
tuations. As in the case of systems in ID, a longer range 
coupling in the form of a power law can maintain long- 
range crystalline order for T > if the decay exponent 
does not exceed a^ D = 3.15. 

Allowing thermally induced fluctuations perpendicular 
to the lattice causes a rapid divergence (i.e. linear) of the 
RMS deviations with the system size, even if the lattice 
geometry is dual-layered in an extended scheme to pro- 
vide local stiffness. The growth of S" ms is asymptotically 
linear even if the coupling between sites is long-ranged, 
decaying as a power law. Hence, long-ranged London in- 
teractions would not be enough by themselves to preserve 
positional order in a covalently bonded locally rigid two 
dimensional lattice. 
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